Create regional pools using BirdLife

source('../env.R')

Generate Worldwide Distributions for all Columbidae

Read in the Birdlife Distribution

Birdlife dataset “Bird Species Distribution Maps of the World.” (2020).

A download of the birdlife distributions can be requested from here: https://datazone.birdlife.org/species/requestdis

BIRDLIFE_DISTRIBUTIONS = '/Users/james/Dropbox/PhD/BirdLife/Distribution/SppDataRequest_columbidae/SppDataRequest.shp'
birdlife = st_read(BIRDLIFE_DISTRIBUTIONS)

Reading layer `SppDataRequest' from data source `/Users/james/Dropbox/PhD/BirdLife/Distribution/SppDataRequest_columbidae/SppDataRequest.shp' using driver `ESRI Shapefile'
Simple feature collection with 529 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -55.72332 xmax: 180 ymax: 70.13647
Geodetic CRS:  WGS 84

head(birdlife)

Simple feature collection with 6 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -46.62738 xmax: 180 ymax: 65.45892
Geodetic CRS:  WGS 84
     SISID            SCI_NAME PRESENCE ORIGIN SEASONAL
1 22690059   Raphus cucullatus        5      1        1
2 22690062 Pezophaps solitaria        5      1        1
3 22690066       Columba livia        1      3        1
4 22690066       Columba livia        1      1        1
5 22690074   Columba rupestris        1      1        3
6 22690074   Columba rupestris        1      1        1
                                                                                                                                          SOURCE
1 Cheke, 1987; Hachisuka, 1953; Livezey, 1993; Owadally, 1979; Strickland and Melville, 1848; Temple, 1977 ;Roberts and Solow, 2004; Cheke, 2006
2                                                                                                     Cheke, 1987; Cowles, 1987; Hachisuka, 1953
3                                            Cramp, 1997; Rasmussen, 2004; Flint, 1984; Ding, 2001; Porter, 1996; Urban, 1986; Gibbs et al, 2001
4                                Cramp, 1997; Rasmussen, 2004; Flint, 1984; Ding, 2001; Porter, 1996; Urban, 1986; Gibbs et al, 2001; eBird 2019
5                                                                                                                          Grimmett et al., 2008
6                                       Gibbs et al., 2001; Rasmussen and Anderton, 2005; Spierenburg, 2005; Grimmett et al., 2008; Brazil, 2009
                                                                        COMPILER DATA_SENS SENS_COMM DIST_COMM TAX_COMM GENERALISD
1                                         Philip Martin (BirdLife International)         0      <NA>      <NA>     <NA>          0
2                                         Philip Martin (BirdLife International)         0      <NA>      <NA>     <NA>          0
3                                                         BirdLife International         0      <NA>      <NA>     <NA>          0
4 Lynda Donaldson (BirdLife International); Donal Smith (BirdLife International)         0      <NA>      <NA>     <NA>          0
5                                            Joe Taylor (BirdLife International)         0      <NA>      <NA>     <NA>          0
6                                            Joe Taylor (BirdLife International)         0      <NA>      <NA>     <NA>          0
                                                       CITATION YRCOMPILED YRMODIFIED VERSION                       geometry
1 BirdLife International and Handbook of the Birds of the World       2010       2016  2022.2 MULTIPOLYGON (((57.64288 -2...
2 BirdLife International and Handbook of the Birds of the World       2010       2016  2022.2 MULTIPOLYGON (((63.47229 -1...
3 BirdLife International and Handbook of the Birds of the World       2006       2016  2022.2 MULTIPOLYGON (((5.383484 53...
4 BirdLife International and Handbook of the Birds of the World       2016       2019  2022.2 MULTIPOLYGON (((88.3183 21....
5 BirdLife International and Handbook of the Birds of the World       2010       2016  2022.2 MULTIPOLYGON (((75.19348 36...
6 BirdLife International and Handbook of the Birds of the World       2010       2016  2022.2 MULTIPOLYGON (((120.8933 36...

Read in the Birdlife Taxonomy

Load the birdlife taxonomy to find Columbidae species.

The taxonomy can be downloaded here: https://datazone.birdlife.org/species/taxonomy

st_crs(birdlife)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

BIRDLIFE_TAXONOMY = '/Users/james/Dropbox/PhD/BirdLife/Taxonomy/Handbook of the Birds of the World and BirdLife International Digital Checklist of the Birds of the World_Version_8.xlsx'
birdlifeTaxonomyWb = read_excel(BIRDLIFE_TAXONOMY)

New names:

birdlifeTaxonomyWb

names(birdlifeTaxonomyWb) <- birdlifeTaxonomyWb[2,]
birdlifeTaxonomy <- birdlifeTaxonomyWb[-2,]
birdlifeTaxonomy

Fetch all columbidae species that are not extinct and are not an invalid birdlife species.

birdlifeColumbidae = birdlifeTaxonomy[!is.na(birdlifeTaxonomy$Order) & birdlifeTaxonomy$Order == 'COLUMBIFORMES', c('SISRecID', 'Scientific name', 'Common name', 'Synonyms', '2023 IUCN Red List category')]
names(birdlifeColumbidae) = c('birdlife_id', 'species_name', 'common_name', 'synonyms', 'red_list_2023')
birdlifeColumbidae = birdlifeColumbidae[!(birdlifeColumbidae$red_list_2023 %in% c('EX', 'NR')),]
birdlifeColumbidae

Check to see if we have any missing distributions

birdlifeColumbidae[!(birdlifeColumbidae$species_name %in% birdlife$SCI_NAME),]

Store Birdlife taxonomy

write_csv(birdlifeColumbidae, filename(BIRDLIFE_WORKING_OUTPUT_DIR, 'taxonomy.csv'))
birdlifeColumbidae = read_csv(filename(BIRDLIFE_WORKING_OUTPUT_DIR, 'taxonomy.csv'))

Rows: 353 Columns: 5── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (4): species_name, common_name, synonyms, red_list_2023
dbl (1): birdlife_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(birdlifeColumbidae)

Filter Birdlife Distributions to only Columbidae

Do we have any distributions not in our taxonomy, are are not extinct?

birdlife[!(birdlife$SCI_NAME %in% birdlifeColumbidae$species_name) & birdlife$PRESENCE != 5,]

Simple feature collection with 0 features and 16 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
 [1] SISID      SCI_NAME   PRESENCE   ORIGIN     SEASONAL   SOURCE     COMPILER   DATA_SENS  SENS_COMM  DIST_COMM  TAX_COMM   GENERALISD CITATION   YRCOMPILED
[15] YRMODIFIED VERSION    geometry  
<0 rows> (or 0-length row.names)

COLUMBIDAE_WORLD_DISTRIBUTIONS = filename(BIRDLIFE_WORKING_OUTPUT_DIR, "columbidae_world.shp")
write_sf(birdlife[birdlife$PRESENCE != 5,], COLUMBIDAE_WORLD_DISTRIBUTIONS)

Warning: GDAL Message 1: Value 154619350 of field SISID of feature 491 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 154619546 of field SISID of feature 492 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 154625310 of field SISID of feature 493 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 154625453 of field SISID of feature 494 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 154625942 of field SISID of feature 495 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 200093329 of field SISID of feature 496 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 200094385 of field SISID of feature 497 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 219582631 of field SISID of feature 498 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 232609896 of field SISID of feature 499 not successfully written. Possibly due to too larger number with respect to field widthWarning: GDAL Message 1: Value 232610475 of field SISID of feature 500 not successfully written. Possibly due to too larger number with respect to field width

columbidaeDistribution = read_sf(COLUMBIDAE_WORLD_DISTRIBUTIONS)

Free up memory by removing birdlife data from workspace

rm(birdlife)
rm(birdlifeTaxonomy)
rm(birdlifeTaxonomy_tmp)

Warning: object 'birdlifeTaxonomy_tmp' not found

rm(birdlifeTaxonomyWb)
rm(birdlifeTaxonomyWb_Sheets)

Warning: object 'birdlifeTaxonomyWb_Sheets' not found

st_crs(columbidaeDistribution)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Create intersection of Columbidae Distributions with Urban Spaces

Read in our initial city vectors, we will join these to our checklists

initial_city_selection = read_sf(filename(mkdir(GEO_WORKING_OUTPUT_DIR, 'cities'), 'initial_selection.shp'))

st_crs(initial_city_selection)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

sf::sf_use_s2(FALSE)

Spherical geometry (s2) switched off

urban_columbidae_distributions = st_join(initial_city_selection, columbidaeDistribution) 

although coordinates are longitude/latitude, st_intersects assumes that they are planar

urban_columbidae_distributions

Simple feature collection with 6595 features and 20 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -158.043 ymin: -38.20685 xmax: 174.9891 ymax: 60.3276
Geodetic CRS:  WGS 84

Check birdlife distribution data

Species with unknown presence?

urban_columbidae_distributions[urban_columbidae_distributions$PRESENCE > 5, c('SCI_NAME', 'city_name', 'city_id')]

Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 116.0291 ymin: -8.698981 xmax: 116.27 ymax: -8.479668
Geodetic CRS:  WGS 84

Possibly Extinct species

urban_columbidae_distributions[urban_columbidae_distributions$PRESENCE == 4, c('SCI_NAME', 'city_name', 'city_id')]

Simple feature collection with 40 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -82.49902 ymin: -25.60243 xmax: 125.6809 ymax: 23.18014
Geodetic CRS:  WGS 84

Map seasonal and origin to string values

urban_columbidae_distributions$seasonal_desc = "Unknown Code"
urban_columbidae_distributions$seasonal_desc[urban_columbidae_distributions$SEASONAL == "1"] = "Resident"
urban_columbidae_distributions$seasonal_desc[urban_columbidae_distributions$SEASONAL == "2"] = "Breeding Season" 
urban_columbidae_distributions$seasonal_desc[urban_columbidae_distributions$SEASONAL == "3"] = "Non-breeding Season" 
urban_columbidae_distributions$seasonal_desc[urban_columbidae_distributions$SEASONAL == "4"] = "Passage" 
urban_columbidae_distributions$seasonal_desc[urban_columbidae_distributions$SEASONAL == "5"] = "Seasonal occurence uncertain"

urban_columbidae_distributions$origin_desc = "Unknown Code"
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "1"] = "Native"
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "2"] = "Reintroduced" 
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "3"] = "Introduced" 
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "4"] = "Vagrant" 
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "5"] = "Origin Uncertain"
urban_columbidae_distributions$origin_desc[urban_columbidae_distributions$ORIGIN == "6"] = "Assisted Colonisation"

urban_columbidae_distributions$presence_desc = "Unknown Code"
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "1"] = "Extant"
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "2"] = "Probably Extant" 
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "3"] = "Possibly Extant" 
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "4"] = "Possibly Extinct" 
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "5"] = "Extinct"
urban_columbidae_distributions$presence_desc[urban_columbidae_distributions$PRESENCE == "6"] = "Presence Uncertain"

cache_urban_columbidae_distributions = filename(BIRDLIFE_WORKING_OUTPUT_DIR, 'urban_distributions_birdlife_only.csv')
write_csv(urban_columbidae_distributions, cache_urban_columbidae_distributions)
urban_columbidae_distributions = read_csv(cache_urban_columbidae_distributions)

Rows: 6595 Columns: 24── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (11): city_name, geometry, SCI_NAME, SOURCE, COMPILER, DIST_COMM, TAX_COMM, CITATION, seasonal_desc, origin_desc, presence_desc
dbl (12): city_id, pop_2015, bu_2015, SISID, PRESENCE, ORIGIN, SEASONAL, DATA_SENS, GENERALISD, YRCOMPILED, YRMODIFIED, VERSION
lgl  (1): SENS_COMM
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Missing North American distribution of Streptopelia decaocto and Columba Livia

Birdlife is missing the recent invasion of Streptopelia decaocto into North America. We supplement the data for Streptopelia decaocto with the range from eBird’s Status and Trends dataset.

This data can be downloaded here: https://science.ebird.org/en/status-and-trends/species/eucdov/downloads

The country boundandaries can be downloaded from the world bank here: https://datacatalog.worldbank.org/search/dataset/0038272/World-Bank-Official-Boundaries

COUNTRY_BOUNDARIES = '/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp'
world_map = st_simplify(st_read(COUNTRY_BOUNDARIES), dTolerance = 0.02)

Reading layer `WB_countries_Admin0_10m' from data source `/Users/james/Dropbox/PhD/WorldBank_countries_Admin0_10m/WB_countries_Admin0_10m.shp' using driver `ESRI Shapefile'
Simple feature collection with 251 features and 52 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -59.47275 xmax: 180 ymax: 83.6341
Geodetic CRS:  WGS 84

Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

Collared Dove

Existing range from Birdlife

sf::sf_use_s2(FALSE)
coldov_vect = st_simplify(columbidaeDistribution[columbidaeDistribution$SCI_NAME == 'Streptopelia decaocto', c('PRESENCE')], dTolerance = 0.02)

Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = coldov_vect, aes(geometry = geometry), fill = "red")

Range from eBird

EBIRD_EUCDOV_RANGE = '/Users/james/Dropbox/PhD/eBird/eucdov_range_2022/eucdov_range_2022.gpkg'
eucdov_range_vect = st_read(EBIRD_EUCDOV_RANGE)

Reading layer `range' from data source `/Users/james/Dropbox/PhD/eBird/eucdov_range_2022/eucdov_range_2022.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -129.1392 ymin: 6.847467 xmax: 128.6997 ymax: 60.44009
Geodetic CRS:  WGS 84

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = eucdov_range_vect["scientific_name"], aes(geometry = geom), fill = "red")

Extract cities from eBird range

sf::sf_use_s2(FALSE)

existing_eucdov_cities = urban_columbidae_distributions$city_id[urban_columbidae_distributions$SCI_NAME == 'Streptopelia decaocto']

eucdov_distribution = st_join(initial_city_selection, eucdov_range_vect) 

although coordinates are longitude/latitude, st_intersects assumes that they are planar

eucdov_distribution[!is.na(eucdov_distribution$species_code) & !(eucdov_distribution$city_id %in% existing_eucdov_cities),]

Simple feature collection with 116 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -123.2852 ymin: 16.69475 xmax: 126.8399 ymax: 51.22359
Geodetic CRS:  WGS 84

missing_eucdov_cities = data.frame(
    species_name = 'Streptopelia decaocto', 
    city_id = eucdov_distribution$city_id[!is.na(eucdov_distribution$species_code) & !(eucdov_distribution$city_id %in% existing_eucdov_cities)],
    city_name= eucdov_distribution$city_name[!is.na(eucdov_distribution$species_code) & !(eucdov_distribution$city_id %in% existing_eucdov_cities)],
    origin = 'Introduced',
    seasonal = 'Resident',
    presence = 'Extant'
)
missing_eucdov_cities

Rock dove

Existing range from Birdlife

sf::sf_use_s2(FALSE)
rocdov_vect = st_simplify(columbidaeDistribution[columbidaeDistribution$SCI_NAME == 'Columba livia', c('PRESENCE')], dTolerance = 0.02)

Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = rocdov_vect, aes(geometry = geometry), fill = "red")

Range from eBird

EBIRD_ROCPIG_RANGE = '/Users/james/Dropbox/PhD/eBird/rocpig_range_2022/rocpig_range_2022.gpkg'
rocpig_range_vect = st_read(EBIRD_ROCPIG_RANGE)

Reading layer `range' from data source `/Users/james/Dropbox/PhD/eBird/rocpig_range_2022/rocpig_range_2022.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -159.9356 ymin: -53.31565 xmax: 178.8437 ymax: 69.31292
Geodetic CRS:  WGS 84

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = rocpig_range_vect["scientific_name"], aes(geometry = geom), fill = "red")

sf::sf_use_s2(FALSE)

existing_rocpig_cities = urban_columbidae_distributions$city_id[urban_columbidae_distributions$SCI_NAME == 'Columba livia']

rocpig_distribution = st_join(initial_city_selection, rocpig_range_vect) 

although coordinates are longitude/latitude, st_intersects assumes that they are planar

rocpig_distribution[!is.na(rocpig_distribution$species_code) & !(rocpig_distribution$city_id %in% existing_rocpig_cities),]

Simple feature collection with 453 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -158.043 ymin: -38.11825 xmax: 141.5939 ymax: 53.66794
Geodetic CRS:  WGS 84

missing_rocpig_cities = data.frame(
    species_name = 'Columba livia', 
    city_id = rocpig_distribution$city_id[!is.na(rocpig_distribution$species_code) & !(rocpig_distribution$city_id %in% existing_rocpig_cities)],
    city_name= rocpig_distribution$city_name[!is.na(rocpig_distribution$species_code) & !(rocpig_distribution$city_id %in% existing_rocpig_cities)],
    origin = 'Introduced',
    seasonal = 'Resident',
    presence = 'Extant'
)
missing_rocpig_cities

Add Spotted dove’s expanded range into China

Existing range from Birdlife

sf::sf_use_s2(FALSE)
spodov_vect = st_simplify(columbidaeDistribution[columbidaeDistribution$SCI_NAME %in% c('Spilopelia chinensis','Spilopelia suratensis','Patagioenas albipennis'), c('SCI_NAME')], dTolerance = 0.02)

Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = spodov_vect, aes(geometry = geometry, fill = SCI_NAME))

Range from eBird

EBIRD_SPODOV_RANGE = '/Users/james/Dropbox/PhD/eBird/spodov_range_2022/spodov_range_2022.gpkg'
spodov_range_vect = st_read(EBIRD_SPODOV_RANGE)

Reading layer `range' from data source `/Users/james/Dropbox/PhD/eBird/spodov_range_2022/spodov_range_2022.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -43.72339 xmax: 180 ymax: 44.8739
Geodetic CRS:  WGS 84

ggplot() +  
  geom_sf(data = world_map, aes(geometry = geometry)) +
  geom_sf(data = spodov_range_vect["scientific_name"], aes(geometry = geom), fill = "red")

sf::sf_use_s2(FALSE)

existing_spodov_cities = urban_columbidae_distributions$city_id[urban_columbidae_distributions$SCI_NAME %in% c('Spilopelia chinensis','Spilopelia suratensis','Patagioenas albipennis')]

spodov_distribution = st_join(initial_city_selection, spodov_range_vect) 

although coordinates are longitude/latitude, st_intersects assumes that they are planar

spodov_distribution[!is.na(spodov_distribution$species_code) & !(spodov_distribution$city_id %in% existing_spodov_cities),]

Simple feature collection with 29 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 70.54817 ymin: 27.96718 xmax: 125.8892 ymax: 43.98709
Geodetic CRS:  WGS 84

missing_spodov_cities = data.frame(
    species_name = 'Spilopelia chinensis', 
    city_id = spodov_distribution$city_id[!is.na(spodov_distribution$species_code) & !(spodov_distribution$city_id %in% existing_spodov_cities)],
    city_name= spodov_distribution$city_name[!is.na(spodov_distribution$species_code) & !(spodov_distribution$city_id %in% existing_spodov_cities)],
    origin = 'Introduced',
    seasonal = 'Resident',
    presence = 'Extant'
)
missing_spodov_cities

Store regional pool result

urban_distribution = rbind(
  data.frame(
    species_name = urban_columbidae_distributions$SCI_NAME, 
    city_id = urban_columbidae_distributions$city_id,
    city_name= urban_columbidae_distributions$city_name,
    origin = urban_columbidae_distributions$origin_desc,
    seasonal = urban_columbidae_distributions$seasonal_desc,
    presence = urban_columbidae_distributions$presence_desc
  ),
 rbind(missing_eucdov_cities, rbind(missing_rocpig_cities, missing_spodov_cities))
)

BIRDLIFE_URBAN_DISTRIBUTION = filename(BIRDLIFE_WORKING_OUTPUT_DIR, 'urban_distributions.csv')
write_csv(urban_distribution, BIRDLIFE_URBAN_DISTRIBUTION)

Examine result

urban_distribution = read_csv(BIRDLIFE_URBAN_DISTRIBUTION)

Rows: 7193 Columns: 6── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): species_name, city_name, origin, seasonal, presence
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Species in Manchester’s regional pool

urban_distribution %>% filter(city_name == 'Manchester') 

Species in Miami’s regional pool

urban_distribution %>% filter(city_name == 'Miami') 

Cities with Common Woodpigeon in regional pool

urban_distribution %>% filter(species_name == 'Columba palumbus')

Cities with Collared Dove in regional pool

urban_distribution %>% filter(species_name == 'Streptopelia decaocto')

Regional pool sizes

urban_distribution %>% group_by(city_id, city_name) %>% summarise(pool_size = n()) %>% arrange(pool_size)

`summarise()` has grouped output by 'city_id'. You can override using the `.groups` argument.

“Bird Species Distribution Maps of the World.” 2020. BirdLife International; Handbook of the Birds of the World. http://datazone.birdlife.org/species/requestdis.